Generation of Higher Harmonics of Monochromatic Electromagnetic Radiation

ABSTRACT

Various methods and systems are provided for generation of higher harmonics of monochromatic electromagnetic radiation. In one embodiment, among others, a system includes a double array of planar structures including a first periodically structured planar structure having a periodic spacing corresponding to a resonant frequency and a second periodically structured planar structure in parallel with the first structure. The system also includes a source configured to illuminate the first structure with a beam of monochromatic electromagnetic radiation at the resonant frequency to produce a higher harmonic of the monochromatic electromagnetic radiation. In another embodiment, a method includes directing a beam of monochromatic electromagnetic radiation at a resonant frequency at a double array of planar structures and generating a higher harmonic of the monochromatic electromagnetic radiation. The double array includes a first periodically structured planar structure and a second periodically structured planar structure in parallel with the first structure.

CROSS REFERENCE TO RELATED APPLICATIONS

This application claims priority to copending U.S. provisional application entitled “GENERATION OF HIGHER HARMONICS OF MONOCHROMATIC ELECTROMAGNETIC RADIATION” having Ser. No. 61/488,971, filed May 23, 2011, the entirety of which is hereby incorporated by reference.

BACKGROUND

The conventional way to generate a second harmonic of light uses a slab of optically non-linear crystal that is illuminated by a laser beam. For example, a beam produced by a diode laser has a wavelength of 980 nm (infrared light). It can be converted into a beam with the wavelength of 490 nm (green light) by transmitting the former through an optically non-linear crystal. The process is used in, e.g., a green-light laser pointer.

BRIEF DESCRIPTION OF THE DRAWINGS

Many aspects of the present disclosure can be better understood with reference to the following drawings. The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views.

FIG. 1A is a cross-sectional view of an example of a double array of dielectric cylinders for generation of higher harmonics of monochromatic electromagnetic radiation in accordance with various embodiments of the present disclosure.

FIG. 1B is a graphical representation illustrating the scattering process exhibited by the cylinders of the double array of FIG. 1A in accordance with various embodiments of the present disclosure.

FIG. 1C is a plot illustrating examples of positions of resonances of the double array of FIG. 1A in accordance with various embodiments of the present disclosure.

FIGS. 2A-2C are plots illustrating the efficiency of second harmonic generation utilizing the double array of FIG. 1A in accordance with various embodiments of the present disclosure.

DETAILED DESCRIPTION

Disclosed herein are various embodiments of methods related to higher harmonic generation of monochromatic electromagnetic radiation. Reference will now be made in detail to the description of the embodiments as illustrated in the drawings, wherein like reference numbers indicate like parts throughout the several views.

Higher harmonics of light may be generated using a slab of optically non-linear material (e.g., a crystal) that is illuminated by a laser beam. For example, as the beam propagates through the slab, some of its radiation may be converted to that of the double frequency. This effect is due to the so called non-linear dielectric susceptibility.

A condition for the process to exist is the phase matching condition. The fundamental field of frequency ω must propagate through the slab of an optically non-linear material in phase with the generated field of the double frequency 2ω, otherwise the amplitude of the second harmonics oscillates and does not grow. The phase mismatch is due to a difference in the propagation speed of light of frequencies ω and 2ω as the refraction index of the slab is typically frequency dependent. There are various techniques to cope with the problem. If the phase matching condition is fulfilled, then the amplitude of the second harmonic grows linearly with the distance traveled through the slab. In a practical set-up, the phase matching condition can only be fulfilled with a certain accuracy. At a certain distance the linear growth stops and no further amplification of the second harmonic is gained. The length (or distance) at which this occurs is an active length specific to the optically non-linear material.

Not all radiation of the incident beam gets converted to a double frequency. After passing through the slab of an active length, the exit beam is a mixture of two beams of frequencies ω and 2ω. Typically a prism placed after the slab is used to separate the single-frequency and double-frequency beams. If P_(ω) is the power per unit area of the incident beam, then the power per unit area of the generated double-frequency beam is P_(2ω)=σ₂P_(ω) where the coefficient 0≦σ₂<1 is the conversion rate. The conversion rate is also specific to an optically non-linear material (assuming the phase matching condition is fulfilled).

A non-linear susceptibility is characterized by constants χ_(n), n=2, 3 . . . for the frequency conversion ω→nω. Its physical significance is that if an electric field E is applied to a material, then the electric field in the material is D=∈E where the dielectric constant ∈ depends on the applied field E:

$\begin{matrix} \begin{matrix} {D = {{ɛ\; E} = {\left( {ɛ_{0} + {\chi_{2}E} + {\chi_{3}E^{2}} + \ldots}\mspace{14mu} \right)E}}} \\ {= {ɛ_{0}\left( {1 + \frac{\chi_{2}E}{ɛ_{0}} + \frac{\chi_{3}E^{2}}{ɛ_{0}} + \ldots}\mspace{14mu} \right)}} \end{matrix} & {{EQN}.\mspace{14mu} 1} \end{matrix}$

where ∈₀ is the linear dielectric constant of the material. The second harmonic generation occurs if the constant χ₂ is not zero for a given material. Although a nonlinear susceptibility χ_(n) is exhibited by many materials, it is typically tiny. As seen from the above equation, the physical parameters that govern optically nonlinear effects are:

$\begin{matrix} {{v_{2} = \frac{\chi_{2}E}{ɛ_{0}}},{v_{3} = {\frac{\chi_{3}E^{2}}{ɛ_{0}}\mspace{14mu} {{etc}.}}}} & {{EQN}.\mspace{14mu} 2} \end{matrix}$

meaning that the only way to amplify the effects is to increase the field amplitude E of the incident beam. In particular, to generate a second harmonic, the parameter ν₂ should be high enough (e.g., ν₂≧0.3). The value achieved in modern devices depends upon the method used to generate the harmonic. For a typical light source, nonlinear effects are highly suppressed (i.e., ν₂<<1 and ν₃ is usually even smaller). Among common light sources, lasers produce beams with the highest field amplitudes. But even laser beams cannot provide for a sufficiently large value of ν₂ in many materials. This problem may be solved by focusing the incident beam to increase power per unit area.

In summary, features to generate higher harmonics include:

-   -   Fulfilling the phase matching condition. This condition is         dictated by the fact that second harmonic generation occurs over         an extend distance traveled by light (the active length), which         is larger than the wavelength of the incident light by about         four to five orders of magnitude.     -   An active length over which a significant conversion rate is         observed varies in a millimeter to centimeter range, while a         typical wavelength of the incident light is on the order of 1000         nm or 10⁻³ mm. So, if L is an active length and λ is the         wavelength of light, then the ratio L/λ is between about 10³ and         10⁴. This parameter is a rough characteristic of the dimensions         of a device to generate higher harmonics. Thus, all such devices         are essentially macroscopic.     -   Focusing of the incident beam can improve the conversion rate         σ=P_(2ω)/P_(ω). A drawback of a tight focusing is that an         optimization of the conversion rate σ becomes a technically         difficult problem.     -   The power conversion rate is typically less than 50%.         Exceptional materials exist for which conversion rates above 80%         have been reported.

Nanophotonics may also be utilized for higher harmonic generation of monochromatic electromagnetic radiation. Periodically structured films with a thickness substantially less than a wavelength of the incident radiation and periods slightly less than the wavelength have been known to exhibit an anomalous transmission and reflection properties, called Wood anomalies. Examples of such periodic structures are a film with a periodic pattern of holes or a film with periodically occurring slits (such structures are also called gratings). Such gratings have been examined in the centimeter and millimeter wave length range. Interest in studying such structures in visible light has increased in the past decade because of the advancements in manufacturing of such periodic structures at a nanoscale, i.e., when the periods of these structures are comparable with the wavelengths of visible light.

The features of scattering of light on periodically structured films are discussed in “Resonant light scattering and higher harmonic generation by periodic sub wavelength arrays” by S. V. Shabanov, International Journal of Modern Physics B, Volume 23, Issue 27, pp. 5191-5236 (2009), which is hereby incorporated by reference in its entirety. The features of scattering of light on periodically structured films include:

First, when the transmitted or reflected power is viewed as a function of the wavelength (or frequency) of the incident light, then it has a narrow peaks at particular wavelength slightly less than the period of the structure. These peaks are called scattering resonances. Their physical origin may be due to periodic geometry of the structure and to physical properties of structure material. Each resonance is characterized by its position, i.e., the wavelength λ_(r) at which the peak occurs, called a resonant wavelength (or, accordingly, a resonant frequency ω_(r)), by its resonant width Γ, which is the size of a frequency interval containing ω_(r) in which the scattered power exceeds half of its maximal value at ω_(r).

Second, when a periodic structure is illuminated with light of the resonant frequency ω_(r), there are regions in the structure in which the local field amplitude E_(r) behaves as:

$\begin{matrix} {E_{r} = \frac{E_{i}}{\sqrt{\Gamma}}} & {{EQN}.\mspace{14mu} 3} \end{matrix}$

where E_(i) is the amplitude of the incident radiation. This implies that a narrow scattering resonance corresponds to high amplification of the field amplitude. In other words, structures that exhibit narrow resonances are capable of focusing the incident radiation at specific regions within the structure. If an optically nonlinear material is placed at such regions, optically nonlinear effects can be significantly amplified. This is especially interesting when the incident radiation is a laser beam, as no focusing is required to amplify optically nonlinear effects generated by laser radiation.

Third, the resonant width Γ may be driven to zero by a mechanism based on the concept of bound states in radiation continuum, which was first noted in quantum mechanics by von Neumann and Wigner. The local field amplification and, hence, optically nonlinear effects may be regulated or controlled using the mechanism. By exploiting the similarity between the basic equations of quantum mechanics and those of electromagnetic theory, it can be shown that the phenomenon exists in electrodynamics. The phenomenon admits the following simple description. Suppose there is a periodically structured film which exhibits a scattering resonance at frequency ω_(r) with width Γ. Consider two such structures placed in parallel at a distance 2h. Using the Breit-Wigner resonant scattering theory, one can prove that:

-   -   The double-array structure has two resonances at frequencies         ω₁(h) and ω₂(h) that depend on the distance h and lie in a         neighborhood of ω_(r); the resonances have widths Γ₁(h) and         Γ₂(h), which also depend on the distance h.     -   There always exists a critical distance h_(b) at which either         the width Γ₁ or Γ₂ vanishes, while the other width becomes 2Γ,         i.e.

Γ₁(h _(b))=0,Γ₂(h _(b))=2Γ or Γ₁(h _(b))=2Γ,Γ₂(h _(b))=0  EQN. 4

-   -   The resonances with the vanishing width correspond to a new type         of waveguide modes localized on the structure. In contrast to         previously known waveguide modes (e.g., those in metal cavities         or in defects in photonic crystals, etc.), they occur in the         spectral range of scattered (diffracted) modes of the light.         This explains the term “the bound states in the radiation         continuum.” When the distance h coincides with a critical value         h_(b), the mode cannot be excited by incident radiation. If h         deviates slightly from h_(b), the existence of such a mode in         the structure causes a sharp scattering resonance with a nearly         vanishing width. This feature is precisely due to the fact that         the mode lies in the spectral range of scattering modes.         Conventional waveguide modes do not lie in the spectral range of         diffracted light and, hence, cannot cause resonances with a         nearly vanishing width.     -   The mechanism to achieve a desired amplification of the local         field amplitude is therefore based on the relation

$\begin{matrix} {{E_{r} = \frac{E_{i}}{\sqrt{\Gamma_{n}\left( {h_{b} + {\Delta \; h}} \right)}}},{\omega = {\omega_{n}\left( {h_{b} + {\Delta \; h}} \right)}},{n = 1},2} & {{EQN}.\mspace{14mu} 5} \end{matrix}$

-   -   where Δh is a small deviation of the distance from the critical         one Δh=h−h_(b). As Δh approaches zero, the amplitude E_(r) in         some regions of the structure grows as much as desired. If the         material of the structure has a nonlinear susceptibility, then         in the regions where the field amplification occurs when         reducing the width Γ_(n), the physical parameters governing the         magnitude of optically nonlinear effects, e.g.,

$\begin{matrix} {v_{2} = {\frac{\chi_{2}E_{r}}{ɛ} = \frac{\chi_{2}E_{i}}{ɛ\sqrt{\Gamma_{n}\left( {h_{b} + {\Delta \; h}} \right)}}}} & {{EQN}.\mspace{14mu} 6} \end{matrix}$

-   -   can be made large enough by taking the distance h sufficiently         close to the critical value h_(b).

The smallest distance at which a resonance with the vanishing occurs is about half of the period of the structure. The effect on a specific example of a double array structure was examined in “Electromagnetic bound states in the radiation continuum for periodic double arrays of sub wavelength dielectric cylinders” by R. F. Ndangali and S. V. Shabanov, Journal of Mathematical Physics, Volume 51, Issue 10, 102901 (2010), which is hereby incorporated by reference in its entirety. Since the period is about the wavelength of the incident radiation, the effect can be achieved at the total thickness of the structure being about the wavelength itself. This allows for enhanced optical non-linear effects at nanoscale and, thus, a miniaturization of optically nonlinear devices.

Referring now to FIG. 1A, shown is a cross-sectional view of an example of a double array 100 of sub wavelength dielectric cylinders 103 for generation of higher harmonics of monochromatic electromagnetic radiation in accordance with various embodiments of the present disclosure. The embodiment of FIG. 1A includes two arrays of parallel dielectric cylinders 103. The cylinders 103 are parallel to the y-axis, which is normal to the plane of FIG. 1A, and indicated by solid circles. The two arrays are separated by a distance 106 of 2h. All distances are measured in the units of the structure period 109 (i.e., the distance between the centers of two neighboring cylinders along the x-axis). The period has been normalized in the example of FIG. 1A. An incident plane wave is characterized by the wave vector k, and the electric field is parallel to the cylinders (the TM polarization). The magnitude of the wave vector k=|k| is related to the frequency as k=ω/c where c is the speed of light. So, the quantity k will be used in place of the frequency ω in the following discussion. The radius R of the cylinders is small as compared to the period, R<<1, and also satisfies the condition Rk<<1, i.e., the wave length λ is about the period of the structure and much larger than the cylinder radius because k=2π/λ so that Rk<<1 implies that R<<λ. For example, with R=0.08D, where D is the structure period (i.e., the distance between two neighboring cylinders in a single array) and λ about D, then R/λ=0.08. The direction of the incident wave is determined by the value of the x-component of the wave vector k_(x), so that every wave in the system is defined by the pair (k, k_(x)). The normal incidence corresponds to k_(x)=0.

If D is the period of the structure, R is the radius of the cylinders, and 2h is the distance between the two arrays, then the conditions that have to be satisfied are: R/D<<1, R/h<<1, and the wavelength λ should be of order D or R/λ<<1 (an example is described in the previous paragraph). The length of cylinders is roughly determined by the incident (laser) beam diameter. The number of cylinders is roughly the beam diameter divided by the structure period. Thus, the geometrical parameters of the system are determined by the incident beam wavelength.

FIG. 1B is a graphical representation of the scattering process exhibited by cylinders 103 of the double array 100. With a normal incidence (k_(x)=0) of the wave, a wave of frequency ω=k/c has two scattering channels 112 indicated by rays with single arrow (one is transmitted and the other is reflected). The generated wave of the double frequency 2ω=2k/c has six scattering channels 115 indicated by double arrows (three transmission channels and three reflection channels).

The system has scattering resonances, as defined above, whose position (frequency) depends of the distance h and the incident angle or k_(x). A resonance occurs at k=k_(r)(h, k_(x)). The solid curve 118 in FIG. 1C illustrates positions of resonances for normal incidence (k_(x)=0) as a function of h, i.e., k=k_(r)(h, 0). The dashed curve of FIG. 1C also corresponds to another set of resonances which will not be used in what follows.

The local field amplification due to a small resonance width, E_(r)=E_(i)/√{square root over (Γ)}, occurs exactly on the cylinders, meaning that if the material of the cylinders has a non-linear susceptibility, the corresponding optical nonlinear effects will be enhanced, provided the width Γ can be made small enough.

The system admits resonances with a vanishing width, which occur at specific distances, as generally described above. The resonances with the vanishing width are indicated by solid dots 121 on the curve k=k_(r)(h, 0). Note that two such states exist in the interval 0<h<1. At least one of the states occurs when the distance between the two arrays does not exceed the structure period (which has been normalized to 1), i.e., when the distance is less than the wavelength of the incident radiation. In the example of FIG. 1C, the first two states occur at about h=0.26 (which corresponds to the distance 2h=0.52) and about h=0.75 (which corresponds to the distance 2h=1.5).

The dielectric response of the material of the cylinders may be characterized by a dielectric constant ∈_(c) and a second-order nonlinear susceptibility χ_(c). The standard nonlinear Maxwell's equations may be solved to find the scattered field of frequency ω and the generated field of frequency 2ω in the so called dipole approximation kR<<1 (which is justified by the geometry of the system). Using the solution, the power P_(2ω) of the generated second harmonic transmitted across a unit area in the direction normal to the structure has been calculated. The conversion efficiency is defined above as σ₂=P_(2ω)/P_(ω), where P_(ω) is the power per unit area of the incident wave. It is a function of the distance h, the radius R, and the material constants ∈ and χ_(c). The efficiency σ₂ has been optimized relative to the distance h. It has been proved that the maximal efficiency occurs when the distance satisfies the condition:

$\begin{matrix} {\left( {h - h_{b}} \right)^{4} = \frac{\chi_{c}^{2}}{8\; \pi^{5}{k_{b}\left( {k_{b}R} \right)}^{6}\left( {ɛ_{c} - 1} \right)}} & {{EQN}.\mspace{14mu} 6} \end{matrix}$

where h_(b) is a distance at which a resonance with the vanishing width occurs (as shown in FIG. 1C), k_(b)=k_(r)(h_(b), 0) is the magnitude of the wave vector k at h=h_(b) (k_(x)=0 for normal incidence), i.e. ω_(b)=ck_(b) is the frequency at which a resonance with the vanishing width OMB's.

The exact expression for the conversion rate σ₂ and its maximal value σ_(2,max) attained at the distance defined by EQN. 6 may be shown as follows. If two identical periodic planar structures that have a scattering resonance at a frequency ω₀ with a width Γ₀ are positioned parallel at some distance 2h, then due to multiple scattering effects the combined system has two resonances at ω_(1,2)=ω_(1,2)(h) with the corresponding widths Γ_(1,2)=Γ_(1,2)(h) which depend on the distance h (and, of course, on ω₀ and Γ₀). For distances large enough to neglect contributions of evanescent fields to the multiple scattering, this dependence can be derived by means of the Breit-Wigner theory. The latter shows, in particular, the existence of localized modes between the two planar structures at some critical values h=h_(b) defined by the conditions Γ₁(h_(b))=0 (whereas Γ₂(h_(b))=2Γ₀) or Γ₂(h_(b))=0 (whereas Γ₁(h_(b))=2Γ₀). For distances small enough so that contributions of evanescent fields in the multiple scattering effects cannot be neglected, a full solution of the scattering problem establishes the existence of bound states in the radiation continuum. For a system of two parallel arrays of periodically positioned sub wavelength dielectric cylinders (e.g., as depicted in FIG. 1A), the existence of bound states in the radiation continuum can be established in numerical studies of the system. A complete classification of bound states as well as their analytic form may be given where bound states exist in higher order diffraction channels, i.e. when more than one decay channel are open, which is quite unusual in resonant scattering theory because by tuning just a single physical parameter h a quasi-stationary state that otherwise decays into several open diffraction channels becomes stationary and decoupled from the radiation continuum.

One of the most interesting features of resonant scattering by planar sub wavelength periodic structures is the existence of “hot” spots within the region occupied by the structure where, for a narrow resonance, the field amplitude E_(r)˜E_(i)/√{square root over (Γ₀)} is magnified as compared to the incident field amplitude E_(i). In particular, for a single periodic array of sub wavelength dielectric cylinders, the field amplification occurs in the cylinders and can be explained by a constructive interference of scattered waves from each cylinder that occurs at the resonant frequency ω₀. If the material of cylinders possesses a nonlinear dielectric susceptibility, then an amplification of optical nonlinear effects are possible.

As noted, in a double-array structure the resonance width can be controlled by varying the distance h in a vicinity of a critical value h_(b) at which a bound state is formed. Some of the regions where the local field amplification E_(r)˜E_(i)/√{square root over (Γ(h)))} occurs in the double array of dielectric cylinders coincide with positions of the cylinders. If χ_(c) is a second order nonlinear dielectric susceptibility, then the quantity χ_(c)E_(r) (that determines the second harmonic generation) can be made large despite the smallness of χ_(c)E_(i) by taking h close enough to its critical value h_(b). This qualitative assessment should be taken with a precaution. A rigorous analysis of the scattering problem by means of the formalism of Siegert states (appropriately extended to periodic structures) shows that no divergence of a local field occurs as h→h_(b). The findings can be summarized as follows.

The analysis is based on the sub wavelength approximation, meaning that the incident wave length is larger than the radius R of cylinders. If k is the magnitude of the wave vector, then the theory has three small parameters:

$\begin{matrix} {{{\delta_{0}(k)} = {{\frac{1}{4}({kR})^{2}\left( {ɛ_{c} - 1} \right)}1}},{\chi_{c}1},{{{\Delta \; h}} = {{{h - h_{b}}}1}}} & {{EQN}.\mspace{14mu} 7} \end{matrix}$

where all the distances are measured in units of the structure period (i.e., R<1), δ₀(k) is the scattering phase for a single cylinder, ∈_(c) is a linear dielectric susceptibility, and the amplitude of the incident wave is set to one, E_(i)=1. With this choice of units, all three parameters are dimensionless. Nonlinear Maxwell's equations are solved by taking into account two nonlinear effects: a second harmonic generation in the leading order in χ_(c) and the fundamental harmonic generation by mixing the second and fundamental harmonics in the leading order in χ_(c) (the second order effect in χ_(c)). The latter ensures the energy flux conservation. In particular for the normal incidence, the ratio of the flux of the second harmonic along the normal direction to the incident flux has been found to be:

σ₂ =Cχ _(c) ² E _(r) ⁴  EQN. 8

where E_(r)=E_(r)(χ_(c),Δh,δ₀) is the field on the cylinders and C=C(δ₀,Δh) is some function. Due to a resonant nature of the scattering process, the function E_(r) is non-analytic in the vicinity of zero values of their arguments so that that the conventional perturbative approach (based on a power series expansion) becomes inapplicable to analyze the conversion rate σ₂. A non-perturbative approach may be used to circumvent the problem. The generated flux of the second harmonic attains its maximal value when the small parameters satisfy the condition:

(Δh)⁴δ₀ ³(k _(b))=μχ_(c) ²  EQN. 9

where μ<<1 is a numerical constant, and k_(b) is the magnitude of the wave vector of the bound state that occurs at h=h_(b). Under this condition, σ₂ becomes analytic in the scattering phase δ₀ so that in the leading order,

σ_(2,max)≈4πk _(b) ⁻¹δ₀(k _(b)).  EQN. 10

An interesting feature to note is the independence of the conversion efficiency of the nonlinear susceptibility χ_(c) (in the leading order in the scattering phase δ₀). In other words, given a nonlinear susceptibility χ_(c), by fine tuning the distance between the arrays it is possible to reach the maximal value determined by the scattering phase at the wave length of a bound state. The largest value of k_(b) for the system considered occurs just below the first diffraction threshold (the wavelength is slightly less than the structure period), i.e. k_(b)≈2π. Taking, for example, R=0.15 and ∈_(c)=2 (so that δ₀(2π) 0.22), the conversion rate reads σ_(2,max)≈0.44, that is, about 44% of the incident flux is converted into the second harmonic flux, which is comparable with the conversion rate achieved in slabs (crystals) of optical nonlinear materials.

The following features of the proposed mechanism to generate higher harmonics are noteworthy. First, the need to fulfill the phase matching condition, which is necessary in optically nonlinear crystals, has been eliminated. The reason is that the second harmonic is generated in regions (cylinders) of dimensions much smaller than the wave length and, hence, the phase mismatch in propagation of the first and second harmonics due to the difference in the corresponding refraction indices does not even occur. Second, the process does not require focusing of the incident beam. Instead, if geometrical parameters of the system are set to maximize the conversion rate, the focusing occurs within the structure automatically and a conversion rate of about 40% is achieved, even though χ_(c)E_(i)<<1 for the incident radiation. Third, the maximum value of the conversion rate depends weakly on the nonlinear susceptibility. This provides a possibility for frequency conversion of lower power light beams. Fourth, an active length at which the conversion rate is maximal is 2h_(b), whose smallest value for the system studied is roughly a half of the wave length of the incident light, while an effective conversion in a slab of an optical nonlinear material requires a length varying between a few millimeters to a few centimeters. This means that the conversion can effectively be done at nanoscales for visible light. These features may be attributed to the existence of resonances with the vanishing width or bound states in the radiation continuum for the scattering system studied, and would not be possible otherwise.

Referring back to FIG. 1A, scattering theory for the system is examined. Consider a system consisting of an infinite double array of parallel, periodically positioned cylinders 103. The cylinders are made of a nonlinear dielectric material with a linear dielectric constant ∈_(c)>1, and a second order susceptibility χ_(c)<<1. For example, cylinder material may include, but is not limited to, lithium niobate or other appropriate dielectric material. The coordinate system is set so that the cylinders are parallel to the y-axis, the structure is periodic along the x-axis, and the z-axis is normal to the structure. The unit of length is taken to be the array period 109, and the distance between the two arrays 106 relative to the period is 2h. The three coordinate axes x, y, and z, are oriented by the unit vectors {e₁, e₂, e₃}. In the case when the structure is illuminated by a plane wave with electric field parallel to the cylinders (TM polarization), Maxwell's equations are reduced to the scalar wave equation,

$\begin{matrix} {{\frac{1}{c^{2}}{\partial_{t}^{2}\left( {{ɛ\; E} + {\frac{\chi}{4\; \pi}E^{2}}} \right)}} = {\Delta \; E}} & {{EQN}.\mspace{14mu} 11} \end{matrix}$

where the dielectric function ∈ has a constant value ∈_(c)>1 on the dielectric cylinders, and equals 1 elsewhere. Similarly, the second order susceptibility χ takes a constant value χ_(c)<<1 on the cylinders and 0 elsewhere. Due to the translation invariance of the scattering structure in the γ-direction and the TM polarization, a solution to EQN. 11 is a function of x and z alone.

It is customary to assume that the solution E is analytic in χ_(c) so that higher harmonic effects may be studied through a power series expansion,

E=2Re{E ₁ e ^(−iwt)+χ_(c) E ₂ e ^(−2iwt)+χ_(c) ²(E _(3,1) e ^(−iwt) +E _(3,3) e ^(−3iwt))+ . . . }  EQN. 12

where E₁ is the amplitude of the fundamental harmonics in the zero order of χ_(c), E₂ is the amplitude of the second harmonics in the first order of χ_(c), and so on. Such a series is obtained by perturbation theory in χ_(c) which entails first solving EQN. 11 in the linear case (χ_(c)=0) to produce a solution E_(L). The general solution E to the nonlinear wave equation is then sought in the form E=E_(L)+E_(NL), where E_(NL) is the correction due to nonlinear effects. If Ĝ is the Green's function of the operator

$\frac{\text{?}}{c^{2}}{\partial_{t}^{2}{- \Delta}}$ ?indicates text missing or illegible when filed  

with appropriate (scattering) boundary conditions, then,

$\begin{matrix} {E_{NL} = {{- \frac{\chi_{c}}{4\pi \; c^{2}}}{\hat{G}\left\lbrack {\eta {\partial_{t}^{2}\left( {E_{L} + E_{NL}} \right)^{2}}} \right\rbrack}}} & {{EQN}.\mspace{14mu} 13} \end{matrix}$

where the second order susceptibility χ has been written in the form χ=χ_(c)η to isolate the perturbation parameter χ_(c). The function η is simply the indicator function of the region occupied by the infinite double array, i.e., its value is 1 on any of the cylinders, and 0 elsewhere. A power series expansion for the correction E_(NL) and, hence, for the full solution E is then obtained by the method of successive approximations.

According to scattering theory, the kernel of the Green's function Ĝ will be meromorphic in

$k^{2} = {\frac{\omega^{2}}{c^{2}}.}$

If the kernel has a real pole k²=k_(b) ²>0, it is not summable, and therefore the successive approximations produce a diverging series, thus indicating a non-analytic behavior of the solution in χ_(c). As is clarified shortly, real poles correspond to bound states in the radiation continuum so that the goal is to investigate the nonlinear wave equation (EQN. 11) in a small neighborhood real poles of the Green's function Ĝ, which amounts to finding a non-analytical solution in χ_(c). This is a crucial difference between conventional treatments of optically nonlinear effects and the present study from both the physical and mathematical points of view. The conclusions summarized in Introduction stem directly from non-analyticity of the solution of the non-linear scattering problem in a spectral region that contains bound states in the radiation continuum.

In general, the problem is posed as a scattering problem. An incident radiation,

E _(in)(r,t)=2 cos(k·r−wt),k=k _(x) e ₁ +k ₂ e 3 ,ck=w,  EQN. 14

is scattered by the double array of dielectric cylinders. The general solution to EQN. 11 should then be of the form,

$\begin{matrix} {{E\left( {r,t} \right)} = {\sum\limits_{l = {- \infty}}^{\infty}{{E_{l}(r)}^{{- }\; l\; \omega \; t}}}} & {{EQN}.\mspace{14mu} 15} \end{matrix}$

where E₀≡0, and for all l, E_(l) is the complex conjugate of E_(l) (as E is real). Therefore it is sufficient to determine only E_(l), l≧1. The amplitudes E_(l) satisfy the Bloch's periodicity condition derived by requiring that the full solution E to the wave equation satisfies the same periodicity condition as the incident wave E_(in), namely,

$\begin{matrix} {{E_{i\; n}\left( {{r + e_{1}},t} \right)} = {E_{i\; n}\left( {r,{t - \frac{k_{x}}{\omega}}} \right)}} & {{EQN}.\mspace{14mu} 16} \end{matrix}$

It then follows that,

E _(l)(r+e ₁)=e ^(ilk) ^(x) E _(l)(r)  EQN. 17

This is the Bloch's condition for the amplitude E_(l). By EQN. 11, the amplitudes of different harmonics satisfy the equations,

$\begin{matrix} {{{{\Delta \; E_{l}} + {l^{2}k^{2}ɛ\; E_{l}}} = {{- {vl}^{2}}{k^{2}\left( {ɛ - 1} \right)}{\sum\limits_{p}{E_{p}E_{{l - p},}}}}}{v = \frac{\chi_{c}}{4\; {\pi \left( {ɛ_{c} - 1} \right)}}}} & {{EQN}.\mspace{14mu} 18} \end{matrix}$

For ease of notation, the parameter ν is often used in lieu of χ_(c). Whenever this is the case, it should be kept in mind that ν is proportional to χ_(c), and is therefore very small.

The scattering theory requires that for 1≠±1, the partial waves E_(l)e^(−i1wt) be outgoing in the spatial infinity (|z|→∞). The fundamental waves E_(±1)e^(∓iwt) are a superposition of an incident plane wave e^(±i(k·r−ωt)) and a scattered wave which is outgoing at the spatial infinity. In all, the above boundary conditions lead to a system of Lippmann-Schwinger integral equations for the amplitudes E_(l):

$\begin{matrix} \left\{ \begin{matrix} {E_{1} = {{{\hat{H}\left( k^{2} \right)}\left\lbrack {E_{1} + {v{\sum\limits_{p}{E_{p}E_{1 - p}}}}} \right\rbrack} + ^{\; {k \cdot r}}}} & \; \\ {{E_{l} = {{\hat{H}\left( ({lk})^{2} \right)}\left\lbrack {E_{l} + {v{\sum\limits_{p}{E_{p}E_{l - p}}}}} \right\rbrack}},} & {l \geq 2} \end{matrix} \right. & {{EQN}.\mspace{14mu} 19} \end{matrix}$

and E_(−l)=Ē_(l) for l≦−1, where Ĥ(q²) is the integral operator defined by the relation:

$\begin{matrix} {{{{\hat{H}\left( q^{2} \right)}\lbrack\psi\rbrack}(r)} = {\frac{q^{2}}{4\; \pi}{\int{\left( {{ɛ\left( r_{0} \right)} - 1} \right){G_{q}\left( {rr_{0}} \right)}{\psi \left( r_{0} \right)}{r_{0}}}}}} & {{EQN}.\mspace{14mu} 20} \end{matrix}$

in which G_(q)(r|r₀) is the Green's function of the Poisson operator, (q²+Δ) G_(q)(r|r₀)=4πδ(r−r₀). For two spatial dimensions, as in the case considered here r=(x,z) and r₀=(x₀,z₀), the Green's function is known to be G_(q)(r|r_(o))=iπH₀(q|r−r₀|) where H₀ is the zero order Hankel function of the first kind.

When ν=0, the amplitudes of all higher harmonics (l≧2) vanish. Therefore it is natural to assume that |E₁|>>|E₂|>>|E₃|>> . . . for a small ν. Note that this does not generally imply that the solution, as a function of ν, is analytic at ν=0. Under this assumption, the solution to the system (EQN. 19) can be approximated by keeping only the leading terms in each of the series involved. In particular, the first equation in EQN. 19 is reduced to

E ₁ ≈e ^(ik·r) +Ĥ(k ²)[E ₁]+2νĤ(k ²)[Ē ₁ E ₂]  EQN. 21

while the second equation becomes

E ₂ ≈Ĥ((2k)²)[E ₂ ]+νĤ((2k)²)[E ₁ ²].  EQN. 22

It then follows that a first order approximation to the solution of the nonlinear wave equation (EQN. 11) may be found by solving the system formed by the EQNS. 21 and 22. To facilitate the subsequent analysis, the system is rewritten as

$\begin{matrix} \left\{ \begin{matrix} {{\left\lbrack {1 - {\hat{H}\left( k^{2} \right)}} \right\rbrack \left\lbrack E_{1} \right\rbrack} = {^{\; {k \cdot r}} + {2\; v{{\hat{H}\left( k^{2} \right)}\left\lbrack {{\overset{\_}{E}}_{1}E_{2}} \right\rbrack}}}} \\ {{\left\lbrack {1 - {\hat{H}\left( \left( {2\; k} \right)^{2} \right)}} \right\rbrack \left\lbrack E_{2} \right\rbrack} = {v\; {{\hat{H}\left( \left( {2\; k} \right)^{2} \right)}\left\lbrack E_{1}^{2} \right\rbrack}}} \end{matrix} \right. & {{EQN}.\mspace{14mu} 23} \end{matrix}$

The solution of the first of EQN. 23 involves inverting the operator 1−Ĥ(k²), and therefore necessitates a study of the poles of the resolvent [1−Ĥ(k²)]¹. Such poles are generalized eigenvalues to the generalized eigenvalue problem,

Ĥ(k ²)[E]=E  EQN. 24

for fixed k_(x). The corresponding eigenfunctions E=E_(s) are referred to as Siegert states. In contrast to Siegert states in quantum scattering theory, electromagnetic Siegert states satisfy the generalized eigenvalue problem (EQN. 24) in which the operator is a nonlinear function of the spectral parameter k². In general, for fixed h, a pole to the resolvent [1−Ĥ(k²)]⁻¹ will have the form k_(r) ²(h)−iΓ(h). Scattered modes satisfy the condition k>k_(x) (the radiation continuum). Hence, if k_(r)(h)>k_(x), then, according to scattering theory, such a pole is a resonance pole. In the case of the linear wave equation, i.e., χ_(c)=ν=0, the scattered flux peaks at k=k_(r)(h) indicating the resonance position, whereas the imaginary part of the pole defines the corresponding resonance width (or a spectral width of the scattered flux peak; a small Γ corresponds to a narrow peak). It appears that there is a discrete set of values of the distance h=h_(b) at which the width of some of resonances vanishes, Γ(h_(b))=0, that is, the resolvent [1−Ĥ(k²)]⁻¹ has a real pole at k_(b)=k_(r)(h_(b))>k_(x). The corresponding eigenfunctions (Siegert states) satisfying EQN. 24 are known as bound states in the radiation continuum. It is the presence of such states that makes the scattering problem defined by the system (EQN. 24) impossible to analyze for real k near k_(r)(h_(b)) by expanding the solution into a power series in χ_(c). It should be noted that there exists another class of bound states in this system for which k²<k_(x) ². These are bound states below the radiation continuum. Such states are not relevant for the present study, and henceforth, bound states are understood as bound states in the radiation continuum.

The objective of what follows is to show that if the parameters (h,k) of the system are chosen appropriately in the vicinity of a critical pair (h_(b),k_(b)) at which a bound state in the radiation continuum forms, a significant portion of the energy flux of the fundamental harmonics is transferred to second harmonics.

Note that as far as the second of EQN. 23 is concerned, it can generally be ensured that if k_(b) ² is an eigenvalue of Ĥ for a bound state, then (2k_(b))² is not an eigenvalue. Consequently, the operator [1−Ĥ((2k)²)]⁻¹ is regular in a neighborhood of k_(b) ². On the other hand, the eigenvalues k_(b) ² are isolated. Therefore for k close to, but not equal to, k_(b), the operator 1−Ĥ(k²) is invertible. It then follows that E₁ satisfies the nonlinear integral equation,

E ₁=(1−Ĥ(k ²))⁻¹ [e ^(ik·r)+2ν² Ĥ(k ²)[Ē ₁(1−Ĥ((2k)²))⁻¹ [Ĥ((2k)²)[E ₁ ²]]]]  EQN. 25

where, in accord with the notation introduced in EQN. 20, the function on which an operator acts is placed in the square brackets following the operator. The operator ν²[1−Ĥ((2k²)]⁻¹ that determines the “nonlinear” part of EQN. 25 cannot be viewed as “small” when k²→k_(b) ² no matter how small ν²≠χ_(c) ² is, which precludes the use of a power series representation of the solution in ν². The behavior of the solution of the integral equation (EQN. 25) will be analyzed for k² near k_(b) ² in the limit of sub wavelength dielectric cylinders. In this limit, the action of the operator Ĥ will be proved to be determined by the action of a 2×2 matrix so that EQN. 25 becomes a system of two quadratic equations. The analysis of this system is considerably simplified by establishing first a particular property of the field ratio η defined as

$\begin{matrix} {{\eta \left( {x,z} \right)} = \frac{E_{1}\left( {x,{- z}} \right)}{E_{1}\left( {x,{+ z}} \right)}} & {{EQN}.\mspace{14mu} 26} \end{matrix}$

As already noted, bound states or, more generally, Siegert states have a specific parity relative to the reflection (x, z)→(x,−z), i.e. they are either odd or even functions of z. This is an immediate consequence of the commutativity of the operator Ĥ and the parity operator {circumflex over (P)} defined by Ĥ[E](x, z)=E(x,−z). Near a pole, the meromorphic expansion of the resolvent [1−Ĥ(k²)]⁻¹ has the form

$\begin{matrix} {E_{1} = {{\frac{\; {C(h)}}{k^{2} - {k_{r}^{2}(h)} + {\; {\Gamma (h)}}}E_{s}} + {O(1)}}} & {{EQN}.\mspace{14mu} 27} \end{matrix}$

where C(h) is some constant depending on h, and E_(s) is an appropriately normalized Siegert state. Consider the curve of resonances

: k=k_(r)(h) in the (h,k)-plane. Along

,

$\begin{matrix} {{E_{1}\left( {x,z} \right)} = {{\frac{C(h)}{\Gamma (h)}{E_{s}\left( {x,z} \right)}} + {O(1)}}} & {{EQN}.\mspace{14mu} 28} \end{matrix}$

Now, as h→h_(b), the width Γ(h) goes to 0, and the Siegert state E_(s) becomes a bound state E_(b) in the radiation continuum. It then follows from EQN. 28 that if C(h) does not go to zero faster than Γ(h) as h→h_(b), i.e. the pole still gives the leading contribution to E₁ in this limit, then

$\begin{matrix} {\left. {\eta \left( {x,z} \right)}\rightarrow\frac{E_{b}\left( {x,{- z}} \right)}{E_{b}\left( {x,{+ z}} \right)} \right. = {\pm 1}} & {{EQN}.\mspace{14mu} 29} \end{matrix}$

depending on whether the bound state E_(b) is even or odd in z. In the case of the linear wave equation (ν=0), the constant C(h) is shown to be proportional to √{square root over (Γ(h))}. Therefore, for a small ν, the assumption that C(h) does not go to zero faster than Γ(h) as h→h_(b) is justified.

The following approach is adopted to solve EQN. 25 near a bound state. First, the curve of resonances

in the (h,k)-plane is found. Next, the equation is solved in terms of the ratio η with the pair (k,h) being on the curve

. The principal part of the amplitude E₁ relative to Δh=h−h_(b) is then evaluated near a critical point (h_(b),k_(b)) on

by taking η to its limit value. This approach reveals a non-analytic dependence of the amplitude E₁ on the small parameters Δh and χ_(c) of the system, which is crucial for the subsequent analysis of the second harmonic generation.

The sub wavelength approximation may be examined as follows. The action of the integral operator Ĥ(q²) in EQN. 25 may be approximated in the limit of sub wavelength dielectric cylinders. The approximation is defined by a small parameter

$\begin{matrix} {{\delta_{0}(q)} = {{\frac{({qR})^{2}}{4}\left( {ɛ_{c} - 1} \right)}1}} & {{EQN}.\mspace{14mu} 30} \end{matrix}$

which is the scattering phase of a plane wave with the wavenumber q on a single cylinder of radius R. For sufficiently small R, this approximation is justified. In this approximation, the integral kernel of Ĥ(q²) is defined by (EQN. 20) and has support on the region occupied by cylinders. The condition (EQN. 30) implies that the wavelength is much larger than the radius R, and therefore field variations within each cylinder may be neglected, ψ(x, z)≈ψ(n,±h) where (n,±h) are positions of the cylinders. The integration in Ĥ(q²)[ψ] yields an infinite sum over positions of the cylinders. By Bloch's condition, ψ(n,±h)=e^(inq) ^(z) ψ(0,±h), so that the function Ĥ(q²)[ψ](x,z) is fully determined by only two values ψ(0,±h). In particular,

Ĥ(q ²)[ψ](0,±h)≈αψ(0,±h)+βψ(0,±h)  EQN. 31

where the coefficients α and β are shown to be:

$\begin{matrix} {{\alpha \left( {q,q_{x}} \right)} = {2\; \pi \; \; {\delta_{0}(q)}\left( {{\sum\limits_{m = {- \infty}}^{\infty}\left( {\frac{1}{q_{z,m}} - \frac{1}{2\; \pi \; {\left( {{m} + 1} \right)}}} \right)} + {\frac{}{\pi}{\ln \left( {2\; \pi \; R} \right)}}} \right)}} & {{EQN}.\mspace{14mu} 32} \\ {\mspace{79mu} {{\beta \left( {q,q_{x},h} \right)} = {2\; \pi \; \; {\delta_{0}(q)}{\sum\limits_{m = {- \infty}}^{\infty}\frac{^{2\; \; {hq}_{z,m}}}{q_{z,m}}}}}} & {{EQN}.\mspace{14mu} 33} \end{matrix}$

where q_(z,m)=√{square root over (q²−(q_(x)+2πm)²)} with the convention that if q²<(q_(x)2πm)², then q_(z,m)=i√{square root over ((q_(x)+2πm)²−q²)}. To obtain the energy flux scattered by the structure, the action of the operator Ĥ(q²) on ψ must be determined in the asymptotic region |z|→∞. It is found that for |z|>h+R,

$\begin{matrix} {{{{\hat{H}\left( q^{2} \right)}\lbrack\psi\rbrack}\left( {x,z} \right)} \approx {2\; \pi \; \; {\delta_{0}(q)}{\sum\limits_{m = {- \infty}}^{\infty}{\left( {{{\psi \left( {0,h} \right)}^{{{z - h}}q_{z,m}}} + {{\psi \left( {0,{- h}} \right)}^{{{z + h}}q_{z,m}}}} \right)\frac{^{\; {x{({q_{x} + {2\; \pi \; m}})}}}}{q_{z,m}}}}}} & {{EQN}.\mspace{14mu} 34} \end{matrix}$

Now that the action of the operator Ĥ(q²) has been established in EQNS. 31 and 34, the amplitudes E₁ and E₂ of the fundamental and second harmonics can be determined by solving the system (EQN. 23). As noted earlier, this will be done along a curve

in the (h,k)-plane defined by k=k_(r)(h) where k_(r)(h) is the real part of a pole of [1−Ĥ(k²)]⁻¹, or equivalently, when the incident radiation has the resonant wave number k=k_(r)(h). This curve is obtained by studying zeros of the determinant det[1−Ĥ(k²)]. With this purpose EQN. 24, that defines Siegert states, is written in the sub wavelength approximation according to the rule (EQN. 31), i.e. the action of Ĥ(k²) is taken at the points (0,±h.). This produces the system,

$\begin{matrix} {{{\left\lbrack {1 -} \right\rbrack \begin{pmatrix} E_{b +} \\ E_{b -} \end{pmatrix}} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}},{= \begin{pmatrix} \alpha & \beta \\ \beta & \alpha \end{pmatrix}}} & {{EQN}.\mspace{14mu} 35} \end{matrix}$

where E_(b±)=E_(b)(0,±h) and the functions α=α(k,k_(x)) and β=β(k,k_(x)) have been defined in the previously. In particular, bound states occur at the points (h_(b),k_(b)) at which the determinant det(1−

) vanishes,

$\begin{matrix} {{\det \begin{pmatrix} {1 - \alpha} & {- \beta} \\ {- \beta} & {1 - \alpha} \end{pmatrix}} = {{\left( {1 - \alpha - \beta} \right)\left( {1 - \alpha + \beta} \right)} = 0}} & {{EQN}.\mspace{14mu} 36} \end{matrix}$

It follows from EQN. 35 that the bound states for which 1−α−β=0 are even in z because E_(b+)=E_(b−) in this case. Similarly, the bound states for which 1−α+β=0 are odd in z. More generally, the poles of the resolvent [1−Ĥ(k²)]⁻¹ are complex zeros of det(1−

). They are found by the conventional scattering theory formalism. Specifically, the resonant wave numbers k²=k_(r) ²(h) are obtained by solving the equation Re{1−α±β}=0 for the spectral parameter k². According to the convention adopted in the representation (EQN. 27), the corresponding resonance widths are defined by

$\begin{matrix} {{{\left\lbrack {1 - \mathcal{H}} \right\rbrack \begin{pmatrix} E_{b +} \\ E_{b -} \end{pmatrix}} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}},{\mathcal{H} = \begin{pmatrix} \alpha & \beta \\ \beta & \alpha \end{pmatrix}}} & {{EQN}.\mspace{14mu} 35} \end{matrix}$

where ∂_(k)z denotes the derivative with respect to k². This definition of the width corresponds to the linearization of Re{1−α±β} near k²=k_(r) ²(h) as a function of k² in the pole factor [1−α±β]⁻¹. The curves of resonances k=k_(r)(h)>0 come in pairs. There is a curve connecting the symmetric bound states and another curve connects the odd ones.

In what follows, only the curve connecting symmetric bound states will be considered. The other curve can be treated similarly. With reference to FIG. 10, the first symmetric bound state occurs when the distance 2h is about half the array period, while the skew-symmetric bound states emerge only at larger distances. The solution obtained near the first symmetric bound state corresponds to the smallest possible transverse dimension of the system (roughly a half of the wave length of the incident radiation) at which a significant enhancement of the second harmonic generation can be observed. Thus, from now on the curve of resonances

refers to the curve in the (h,k)-plane defined by the equation Re{1−α−β}=0. To simplify the technicalities, it will be further assumed that only one diffraction channel is open for the fundamental harmonics, i.e. k_(x)<k<2π−k_(x). Note that bound states in the radiation continuum exist even if more diffraction channels are open, and the solution of the nonlinear scattering problem may also be obtained in their vicinity by the approach developed here.

Let k=k_(b)(h) be a solution of Re{1−α−β}=0. By making use of the explicit form of the functions α and β for one open diffraction channel, one infers that along the curve k=k_(b)(h):

$\begin{matrix} \begin{matrix} {{1 - \alpha - \beta} = {i\; {Im}\left\{ {1 - \alpha - \beta} \right\}}} \\ {{= {{- i}\frac{4\; \pi \; {\delta_{0}\left( k_{b} \right)}}{k_{z}}\phi^{2}}},} \end{matrix} & {{EQN}.\mspace{14mu} 38} \\ {{\phi = {\cos \left( {hk}_{z} \right)}},} & \; \\ {k_{z} = \sqrt{k_{b}^{2} - k_{x}^{2}}} & \; \end{matrix}$

Bound states in the radiation continuum occur when the distance h satisfies the equation φ=0, i.e. h√{square root over (k_(r) ²(h)−k_(x) ²)}=(n−1/2)π with n being a positive integer. Its solutions h=h_(b)(n) define the corresponding values of the wave numbers of the bound states, k_(b)(n)=k_(b)(h_(b)(n)). So, the sequence of pairs [h_(b)(n), k_(b)(n))]_(n=1) ^(∝) indicates positions of the bound states on the curve

. In the limit h→h_(b)(n) along

, the function φ has the asymptotic behavior:

φ=(−1)^(n) k _(z,b) Δh+o(Δh),k _(z,b) =k _(z)|_(h=h) _(b(n)) ,Δh=h−h _(b)  EQN. 39

The objective is to determine the dependence of the amplitudes E₁ and E₂ on the parameters Δh and χ_(c) which are both small, and, in particular, to investigate the outgoing flux of the second harmonics as a function of Δh and χ_(c).

To this end, let E_(1±)=E₁(0,±h) be the values of the field E₁ on the axes of the cylinders at (0,±h), and E_(2±)=E₂(0,±h) be the values of the field E₂ on the same cylinders. In the sub wavelength approximation, these values fully determine the scattered field as is shown later and, hence, have to be found first. Applying the rule (EQN. 31) to evaluate the action of the operator Ĥ(q²) in the system (EQN. 23), the first equation of the latter becomes

$\begin{matrix} {{\left\lbrack {1 - \mathcal{H}} \right\rbrack \begin{pmatrix} E_{1 +} \\ E_{1 -} \end{pmatrix}} = {{2\; \upsilon \; {\mathcal{H}\begin{pmatrix} {{\overset{\_}{E}}_{1 +}E_{2 +}} \\ {{\overset{\_}{E}}_{1 -}E_{2 -}} \end{pmatrix}}} + \begin{pmatrix} e^{+ {ihk}_{z}} \\ e^{- {ihk}_{z}} \end{pmatrix}}} & {{EQN}.\mspace{14mu} 40} \end{matrix}$

Similarly, the second of EQN. 23 yields,

$\begin{matrix} {{{\left\lbrack {1 - \mathcal{H}_{2}} \right\rbrack \begin{pmatrix} E_{2 +} \\ E_{2 -} \end{pmatrix}} = {\upsilon \; {\mathcal{H}_{2}\begin{pmatrix} E_{1 +}^{2} \\ E_{1 -}^{2} \end{pmatrix}}}},{\mathcal{H}_{2} = {\mathcal{H}\left( {{2k},{2k_{x}}} \right)}}} & {{EQN}.\mspace{14mu} 41} \end{matrix}$

As stated above, the resolvent [1−Ĥ((2k)²)]⁻¹ is regular in a neighborhood of k_(b) so that EQN. 41 can be solved for E_(2±), which defines the latter as functions of E_(1±). The substitution of this solution into EQN. 40 gives a system of two nonlinear equations for the fields E_(1±). Adding these equations and replacing the field E¹⁻ by its expression E¹⁻=η(0, h)E₁₊ in terms of the field ratio of EQN. 26 yields the following implicit relation between the field E₁₊ and its amplitude |E₁₊|:

$\begin{matrix} {E_{1 +} = {- \frac{\phi}{{\frac{\upsilon^{2}}{\zeta}{E_{1 +}}^{2}} + {\phi^{2}\xi}}}} & {{EQN}.\mspace{14mu} 42} \end{matrix}$

where φ=cos(hk_(z)) and

$\upsilon = \frac{Xc}{4\; {\pi \left( {s_{c} - 1} \right)}}$

are small and, in terms of the field ratio η=η(0, h), the values of ζ and ξ read:

$\begin{matrix} {{\xi = {i\frac{2\; \pi \; {\delta_{0}(k)}}{k_{z}}\left( {1 + \eta} \right)}},} & {{EQN}.\mspace{14mu} 43} \\ {\frac{1}{\zeta} = {\left( {1 + {i\frac{4\; \pi \; {\delta_{0}(k)}}{k_{z}}\phi^{2}}} \right)\left( {a + {b\; \eta^{2}} + {\overset{\_}{\eta}\left( {b + {a\; \eta^{2}}} \right)}} \right)}} & \; \end{matrix}$

with a and b being defined by the relation:

$\begin{matrix} {{\left\lbrack {1 - \mathcal{H}_{2}} \right\rbrack^{- 1}\mathcal{H}_{2}} = \begin{pmatrix} a & b \\ b & a \end{pmatrix}} & {{EQN}.\mspace{14mu} 44} \end{matrix}$

In particular, ζ and ξ are continuous functions of the field ratio η and φ. It can be shown that if ζ_(b) and ξ_(b) are the respective limits of ζ and ξ as h→h_(b) along the curve of resonances

, then these limits are nonzero. It follows then that EQN. 42 for E₁₊ is singular in both η and φ when these parameters are small, i.e., in the limit (η,φ)→(0, 0). Furthermore, there is no way to solve the said equation perturbatively in either of the parameters. A full non-perturbative solution can be obtained using Cardano's method for solving cubic polynomials. Indeed, by taking the modulus squared of both sides of the equation, it is found that,

$\begin{matrix} {{{X^{3} + {2\frac{\phi^{2}}{v^{2}}{Re}\left\{ {\zeta\xi} \right\} X^{2}} + {\frac{\phi^{4}}{v^{4}}{{\zeta\xi}}^{2}X} - {\frac{\phi^{2}}{v^{4}}{\zeta }^{2}}} = 0},{X = {{E_{1 +}}^{2}.}}} & {{EQN}.\mspace{14mu} 45} \end{matrix}$

The solution to this cubic equation may be obtained, where EQN. (24) admits a unique real solution so that there is no ambiguity on the choice of E₁₊. In the vicinity of a point (h_(b),k_(b)) along the resonance curve

, the field E₁₊ is found to behave as,

$\begin{matrix} {{E_{1 +}} = {\frac{{{\Delta \; h}}^{1/3}}{\chi_{c}}{\tau \left( {{\Delta \; h},\chi_{c}} \right)}}} & {{EQN}.\mspace{14mu} 46} \end{matrix}$

Recall that Δh=h−h_(b). An explicit form of the function τ(Δh,χ_(c)) may be derived. It involves combinations of the square and cube roots of functions in Ah and x, and has the property that τ(Δh,χ_(c))→0 as (Δh,χ_(c))→(0, 0) (in the sense of the two-dimensional limit). In the limit h→h_(b), the field ratio η approaches 1 for a symmetric bound state as argued earlier. Therefore it follows from EQN. 41 that E_(2±)νE₁₊ ² because the matrix (EQN. 44) exists at h=h_(b). Since ν≠χ_(c), relation (EQN. 43) leads to the conclusion that:

$\begin{matrix} {\frac{E_{2 \pm}}{E_{1 +}} = {O\left( {{\Delta \; h}}^{1/3} \right)}} & {{EQN}.\mspace{14mu} 47} \end{matrix}$

Thus, the approximation |E₁|>>|E₂| used to truncate the system (EQN. 23) remains valid for h close to the critical value h_(b) despite of non-analyticity of the amplitudes at (Δh,χ_(c))=(0, 0).

The conversion efficiency may be examined using flux analysis. For the nonlinear system considered, even though Poynting's theorem takes a slightly different form as compared to linear Maxwell's equations, the flux conservation for the time averaged Poynting vector holds. The scattered energy flux carried across a closed surface by each of the different harmonics adds up to the incident flux across that surface. Consider a closed surface that consists of four faces,

$L_{\pm} = \left\{ {{\left( {x,z} \right){{- \frac{1}{2}} \leq x \leq \frac{1}{2}}},\left. z\rightarrow{{\pm 1}/2} \right.} \right\}$

and L_(±1/2)={(χ,z)|χ=±1/2}L±1/2={(x, z)|x=±1/2} as depicted in FIG. 1B. The scattered flux of each l^(th)-harmonics across the union of the faces L_(±1/2) vanishes because of the Bloch condition (and so does the incident flux for any k_(x)). Therefore only the flux conservation across the union of the faces L_(±) has to be analyzed. If σ_(l) designates the ratio of the scattered flux carried by the l^(th)-harmonics across the faces L_(±) to the incident flux across the same faces, then Σ_(l≦1)σ_(l)=1. Thus, for l≧2, σ_(l) defines the conversion ratio of fundamental harmonics into the l^(th)-harmonics.

In the perturbation theory used here, only the ratios σ₁ and σ₂ may be evaluated. It can be shown that σ₁+σ₂≦1. Hence, the efficiency of converting the fundamental harmonic into the second harmonic is simply determined by the maximum value of σ₂ as a function of the parameter h at a given value of the nonlinear susceptibility χ_(c).

The ratio σ₂ is defined in terms of the scattering amplitudes of the second harmonic, i.e. by the amplitude of E₂ in the asymptotic region |z|→∞:

$\begin{matrix} \left. {E_{2}(r)}\rightarrow\left\{ \begin{matrix} {{\sum\limits_{m^{{op},{sh}}}{R_{m}^{sh}^{{r} \cdot k_{m,{sh}}^{-}}}},} & \left. z\rightarrow{- \infty} \right. \\ {{\sum\limits_{m^{{op},{sh}}}{T_{m}^{sh}^{\; {r \cdot k_{m,{sh}}^{+}}}}},} & \left. z\rightarrow\infty \right. \end{matrix} \right. \right. & {{EQN}.\mspace{14mu} 48} \end{matrix}$

where k_(m,sh) ^(∓)=(2k_(x)+2πm)e₁±k_(z,rn) ^(sk)e₃ is the wave vector of the second harmonic in the m^(th) open diffraction channel. Recall that the m^(th) channel is open provided (2k)²>(2k_(x)+27πm)² and in this case k_(z,m) ^(sh)=√{square root over ((2k)²−(2k_(x)+2πm)²)}{square root over ((2k)²−(2k_(x)+2πm)²)}, while if the channel is closed, then k_(z,m) ^(sh)=i√{square root over ((2k_(x)+2πm)²−(2k)²)}{square root over ((2k_(x)+2πm)²−(2k)²)}. In the asymptotic region |z|→∞, the field in closed channels decays exponentially and, hence, the energy flux can only be carried in open channels to the spatial infinity. The summation in EQN. 48 is taken only over those values of m for which the corresponding diffraction channel is open for the second harmonic, which is indicated by the superscript “op,sh” in the summation index m^(op,sh). Note that there is more than one open diffraction channel for the second harmonic even though only one diffraction channel is open for the fundamental one. For instance, if the x-component of the wave vector k, i.e. k_(x), is less than

$\frac{\pi}{2},$

there are 3 open diffraction channels for the second harmonic, the channels m=0, m=−1, and m=1. These three directions of the wave vector of the second harmonic propagating in each of the asymptotic regions z−±∞ are depicted in FIG. 1B by double-arrow rays. Thus, in terms of the scattering amplitudes introduced in EQN. 48, the ratio of the second harmonic flux across L_(±) to the incident flux is:

$\begin{matrix} {\sigma_{2} = {\frac{1}{2k_{z}}{\sum\limits_{m^{{op},{sh}}}{{k_{z,m}^{sh}\left( {{R_{m}^{sh}}^{2} + {T_{m}^{sh}}^{2}} \right)}.}}}} & {{EQN}.\mspace{14mu} 49} \end{matrix}$

The scattering amplitudes R_(rn) ^(sh) and T_(rn) ^(sh) may be inferred from EQN. 22 in which the rule (EQN. 34) is applied to calculate the action of the operator Ĥ((2k)²) in the far-field regions |z|→∞:

$\begin{matrix} \left\{ \begin{matrix} {R_{m}^{sh} = {\frac{2\pi \; \; {\delta_{0}\left( {2k} \right)}}{k_{z,m}^{sh}}\left\lbrack {{\left( {E_{2 +} + {vE}_{1 +}^{2}} \right)^{\; {hk}_{z,m}^{sh}}} + {\left( {E_{2 -} + {vE}_{1 -}^{2}} \right)^{{- }\; {hk}_{z,m}^{sh}}}} \right\rbrack}} \\ {T_{m}^{sh} = {{\frac{2{\pi }\; {\delta_{0}\left( {2k} \right)}}{k_{z,m}^{sh}}\left\lbrack {{\left( {E_{2 +} + {vE}_{1 +}^{2}} \right)^{{- }\; {hk}_{z,m}^{sh}}} + {\left( {E_{2 -} + {vE}_{1 -}^{2}} \right)^{\; {hk}_{z,m}^{sh}}}} \right\rbrack}.}} \end{matrix} \right. & {{EQN}.\mspace{14mu} 50} \end{matrix}$

Since η→1 as h→h_(b) along

, the principal part of σ₂ in a vicinity of a bound state along the curve

is obtained by setting E¹⁻=E₁₊ in EQN. 41, solving the latter for E_(2±), and substituting the solution into the expression for σ₂. The result reads:

σ₂ =C _(b)ν² |E ₁₊|⁴  EQN. 51

where C_(b) is a constant obtained by taking all nonsingular factors in the expression of σ₂ to their limit as h→h_(b), which gives

                                   EQN.  52 $C_{b} = \left\lbrack {\frac{\left( {16{{\pi\delta}_{0}(k)}} \right)^{2}}{k_{z}}{{1 + a + b}}^{2}{\sum\limits_{m^{{op},{sh}}}\frac{\cos^{2}\left( {hk}_{z,m}^{sh} \right)}{k_{z,m}^{sh}}}} \right\rbrack_{{({h,k})} = {({h_{b},k_{b}})}}$

for a and b defined in EQN. 44. Using the identity |E₁₊|⁴=|E₁₊|²|E₁₊|² and substituting EQN. 42 into one of the factors |E₁₊|², the conversion ratio σ₂ is expressed as a function of a single real variable,

$\begin{matrix} {{{\sigma_{2}(u)} = {C_{b}^{\prime}\frac{u}{{{u + {\zeta_{b}\xi_{b}}}}^{2}}}},{u = \left( \frac{v{E_{1 +}}}{\phi} \right)^{2}}} & {{EQN}.\mspace{14mu} 53} \end{matrix}$

where C_(b) ^(t)=C_(b)|ζ_(b)|² is a constant, and ξ and ζ in EQN. 43 have been taken at their limits as h→h_(b) to obtain the principal part of σ₂. The function u→σ₂(u) on [0,∞) is found to attain its absolute maximum at u=|ξ_(b)ζ_(b)|. This condition determines the distance 2h between the arrays at which the conversion rate is maximal for given parameters R, ∈_(c) and χ_(c) of the system. Indeed, since u=|E₁₊|² should also satisfy the cubic equation (EQN. 45), the substitution of u=|ξ_(b)ζ_(b)| into the latter yields the condition

$\begin{matrix} {\frac{v^{2}}{\phi^{4}} = {2{\xi_{b}}^{2}\left( {{{\xi_{b}\zeta_{b}}} + {{Re}\left\{ {\xi_{b}\zeta_{b}} \right\}}} \right)}} & {{EQN}.\mspace{14mu} 54} \end{matrix}$

In particular, in the leading order in δ₀(k), the optimal distance 2h between the two arrays is given by the formula,

$\begin{matrix} {\left( {h - h_{b}} \right)^{4} = \frac{\chi_{c}^{2}}{8\pi^{5}{k_{z,b}\left( {k_{b}R} \right)}^{6}\left( {ɛ_{c} - 1} \right)^{5}}} & {{EQN}.\mspace{14mu} 55} \end{matrix}$

where as previously, k_(z,b)=√{square root over (k_(b) ²−k_(x) ²)}.

The maximum value σ_(2,max) of the conversion ratio σ₂ is the sought-for conversion efficiency. Note that σ_(2,max)=σ₂(|ξ_(b)ζ_(b)) is independent of the nonlinear susceptibility χ_(c) because the constants C′_(b), ξ_(b), and ζ_(b) are fully determined by the position of the bound state (k_(b),h_(b)). In other words, if the distance between the arrays is chosen to satisfy the condition (EQN. 55), the conversion efficiency σ_(2,max) is the same for a wide range of values of the nonlinear susceptibility χc. This conclusion follows from two assumptions made in the analysis. First, the sub wavelength approximation should be valid for both the fundamental and second harmonics, i.e. the radius of cylinders should be small enough. Second, the values of h−h_(b) and ν (or x_(c)) must be such that the analysis of the existence and uniqueness of |E₁₊| holds, that is, EQN. 45 should have a unique real solution under the condition (EQN. 55). The geometrical and physical parameters of the studied system can always be chosen to justify these two assumptions as illustrated in FIGS. 2A-2C. A detailed derivation is provided in Appendix A, entitled “The resonant nonlinear scattering theory with bound states in the radiation continuum and the second harmonic generation,” which is hereby incorporated by reference in its entirety.

Referring now to FIGS. 2A-2C, shown are examples of plots illustrating the efficiency of second harmonic generation utilizing the double array 100 of sub wavelength dielectric cylinders 103 of FIG. 1A in accordance with various embodiments of the present disclosure. FIG. 2A shows the maximal efficiency σ_(2,max) as function of R for the case of normal incidence (k_(x)=0). The three curves 203 corresponds to three first critical values h_(b) (i.e., the three first solid dots 121 from the left on the curve 118 in FIG. 1C) and ∈_(c)=1.5. The dashed regions on the curves 203 indicate regions where the dipole approximation kR<<1 used in the analysis was not justified. As one can see, the maximum convergence above 40% occur for values of the radius at which the approximation kR <<1 is well justified. FIG. 2B shows three curves 306 for the maximum efficiency as a function of k_(x) (not normal incidence) for the same three first critical values h_(b) with R=0.15 and ∈_(c)=1.5. As can be seen, the conversion rate of between 35% and 40% can be achieved for a wide range of incident angles. FIG. 2C shows a shaded region of values (ν, h), where ν=χ₂/[4π(∈_(c)6−1)] for which the presented analysis is valid near the second critical value h_(b) and the normal incidence (k, =0). The regions of validity look for other critical values of h look very similar. The found value of h near h_(b) at which the efficiency σ₂ attains its maximal value lies within the shadowed region.

Referring to FIGS. 2A and 2B, shown is the conversion efficiency σ_(2,max) for three first symmetric bound states n=1, 2, 3, as, respectively, a function of the cylinder radius R when k_(x)=0 and of k_(x) when R=0.15. For all curves presented in the panels, ∈_(c)=1.5. The values of σ_(2,max) are evaluated numerically by EQN. 53 where u=|ξ_(b)ζ_(b)ζ. The solid parts of the curves in FIG. 2A correspond to the scattering phase δ₀(k)<0.25 with k=k_(b). Note that the wavelength at which the second harmonic generation is most efficient is the resonant wavelength defined by k=k_(r)(h) where h satisfies the condition (EQN. 55). For a small χ_(c), the scattering phase at the resonant wavelength can well be approximated as δ₀(k)≈∈₀(k_(b)). The condition δ₀(k_(b))<0.25 ensures that the scattering phase for the second harmonic satisfies the inequality δ₀(2k_(b))=4δ₀(k_(b))<1 otherwise the validity of the sub wavelength approximation cannot be justified. The dashed parts of the curves in FIG. 2A correspond to the region where δ₀(k_(b))>0.25. FIGS. 2A and 2B show that the conversion efficiency can be as high as 40% for a wide range of the incident angles and values of the cylinder radius. Such a conversion efficiency is comparable with that achieved in optically nonlinear crystals at a typical beam propagation length (active length) of a few centimeters, whereas here the transverse dimension 2h of the system studied here can be as low as a half of the wavelength, i.e. for an infrared incident radiation, 2h is about a few hundred nanometers. Indeed, as one see in FIG. 1C, the first bound state occurs at h_(b)(1)≈0.259 and k_(b)≈2π which corresponds to the wavelength λ_(b) _(—) b=2π/k_(b)≈1.

The stated conversion efficiency can be fairly well estimated in the leading order of δ₀(k):

$\begin{matrix} {{\sigma_{2,\max} = {{\sigma_{2}\left( {{\zeta_{b}\xi_{b}}} \right)} \approx {8{{\pi\delta}_{0}(k)}{\sum\limits_{m^{{op},{sh}}}\frac{\cos^{2}\left( {hk}_{z,m}^{sh} \right)}{k_{z,m}^{sh}}}}}}}_{{({h,k})} = {({h_{b},k_{b}})}} & {{EQN}.\mspace{14mu} 56} \end{matrix}$

Suppose that only one diffraction channel is open for the incident radiation. Then m=0,±1 in EQN. 56 (three open channels for the second harmonics). Let σ_(2,max) ⁰ denote the term m=0, i.e. σ_(2,max) ⁰ is the second harmonic flux in which the contribution of the channels with m=±1 is omitted, σ_(2,max) ⁰<σ_(2,max). One can infer from EQN. 56 that:

$\begin{matrix} {\sigma_{2,\max}^{0} \approx {\frac{4{{\pi\delta}_{0}\left( k_{b} \right)}}{k_{z,b}}{\cos^{2}\left( {2h_{b}k_{z,b}} \right)}}} & {{EQN}.\mspace{14mu} 57} \end{matrix}$

As the pair (h_(b),k_(b)) at which a symmetric bound state is formed satisfies the equation cos(h_(b)k_(z,b))=0, it follows that cos(2h_(b)k_(z,b))=−1. Hence,

$\begin{matrix} {\sigma_{2,\max}^{0} \approx {\frac{k_{b}^{2}}{k_{z,b}}\pi \; {R^{2}\left( {ɛ_{c} - 1} \right)}}} & {{EQN}.\mspace{14mu} 58} \end{matrix}$

The wavenumbers k_(b) at which the bound states occur lie just below the diffraction threshold 2π−k_(b), i.e., k_(b)≦2π−k_(x) (see FIG. 1C). So that in the case of normal incidence (k_(x)=0), the above estimate becomes,

σ_(2,max) ⁰≈2π(πR ²)(∈_(c)−1)  EQN. 59

with δ₀(2π)<<1. If, for instance, R=0.15 and ∈_(c)=2, then,

σ_(2,max) ⁰≈44%,  EQN. 60

for δ₀(2π)≈0.22.

Apart from the exact numerical results presented in FIGS. 2A and 2B, the main properties of the proposed system to generate higher harmonics can also be understood in the leading order of the dipole (or sub wavelength) approximation. In the approximation, the maximum value of the efficiency is found to be

σ_(2,max) ≈πk _(b) R ²(∈_(c)−1)  EQN. 61

Note that the maximum value of σ₂ appears independent of the non-linear susceptibility. This proves that the conversion efficiency depends weakly on the nonlinear susceptibility, i.e., only through higher orders of the dipole approximation. The only assumption used to prove this fact is that χ_(c)>0 (i.e., some small positive number). This is rather counterintuitive and in a drastic contrast with the conventional approach to the second harmonic generation where σ₂ is always proportional to (χ_(c)E_(i))², where E_(i) is the amplitude of the incident wave. As a point of fact, the efficiency σ₂ in the proposed approach is proved to have a similar structure, namely,

σ₂φ(χ₂ E _(r))²  EQN. 62

where φ=φ(R, h) is some numerical factor that depends on h, R, and ∈_(c), and E_(r) is the electric field on the cylinders (i.e., not the amplitude of the incident wave). This is what makes so much difference with the conventional approach. From the point of view of the resonant scattering theory, it also makes a perfect sense.

First, as noted above, E_(r)=E_(i)/√{square root over (Γ(h))}. So, as h approaches a critical value h_(b), the width Γ(h) tends to zero and the field E_(r) diverges. The system exhibits much needed focusing to amplify non-linear effects, and no special efforts are required to achieve the focusing in full contrast to the conventional case. Second, the physical (or observed) quantity is the flux of scattered light (or the power transmitted across unit area in a specified direction). It remains finite despite the divergence of E_(r). This happens because the function φ(R, h) vanishes as h approaches h_(b) at a faster rate than E_(r) diverges, so that σ₂ remains finite near h=h_(b) and attains its maximum value at h near h_(b) as defined in EQN. 6. Third, the value of h=h_(max) at which the maximum of σ₂ occurs depends on χ_(c), h_(max)=h_(max)(χ_(c)) so that σ₂˜χ_(c) ²/Γ(h_(max)(χ_(c))) is independent of χ_(c) in the leading order of the dipole approximation or, generally, depends weakly on it as the leading order gives a dominant contribution to the conversion efficiency for the specified geometrical parameters of the system.

Furthermore, as illustrated in FIG. 1C, k_(b) are close to 2π (6.05<k_(b)<2π≈6.28). For typical transparent materials ∈≈2. Taking k_(b)=2π, ∈_(c)=2, and R=0.15 (in units of the structure period), for which the dipole approximation is well justified, as noted above the conversion efficiency reads:

σ_(2,max)≈2π(πR ²)(∈_(c)−1)·100%=44%  EQN. 63

Thus, the system exhibits efficiencies consistent with the conventional one, but in contrast to the latter the conversion occurs at distances of order of the wavelength of the incident light. The wavelength is typically of order 1000 nm or 10⁻³ mm versus a typical length of a 10 mm crystal used in the conventional way to generate higher harmonics. There is a distance between two periodic arrays (as described above) at which the second harmonic generation is significantly amplified and reaches or exceeds the level achieved in slabs (crystals) of optically nonlinear materials with a size of about 20 to 30 mm along the beam propagation direction that are currently used in industry. This distance is proved to be of the order of a wavelength of the resonant frequency, thus making the size of the structure along the beam propagation direction about 10000 times smaller than the typical size of crystals.

Compared to conventional use of crystal slabs used to generate the second harmonic, the following advantages have been illustrated. The phase matching condition has been eliminated. The reason is that the second harmonic is generated in regions of the cylinders of dimensions much smaller than the wave length and, hence, the phase mismatch in propagation of the first and second harmonics does not even occur. Focusing of the incident beam is not needed. If geometrical parameters of the system are set to maximize the conversion rate, the focusing occurs within the structure automatically. The maximum value of the conversion rate depends weakly on the nonlinear susceptibility in the sense explained above. So materials with smaller values of the nonlinear susceptibility may be used to achieve the effect. Additionally, this feature provides a possibility for the frequency conversion of lower power light beams, a non-existing option in the conventional approaches. The active length at which the conversion rate is maximal is of the order of the wavelength of the incident beam, reducing it by roughly 10⁴ times. Thus, in contrast to the conventional way to generate higher harmonics, the active length lies in the nanoscale range for visible and infrared light, allowing for miniaturization of devices for higher harmonic generation. The disclosed system allows for a controlled enhancement of local field amplitudes by pairing two nanoscale periodic arrays at specific distances.

Other forms of double periodic structures may also be used to enhance and control nonlinear optical phenomena at the nanoscale. While the effects of a double array of cylinders have been described, other periodic structures, such as periodically structured planar (film-like) structures, when paired at a certain distance can exhibit similar optical properties as those described for the double arrays of cylinders. This distance, at which the enhancement is maximal, depends on the planar structure geometry. Thus, the pairing of such periodic structures at a specific distance can significantly enhance optically nonlinear effects such as, e.g., higher (second) harmonic generation.

Further, the relative position of the double arrays of cylinders may be adjustable. For example, microelectromechanical systems may be used to displace the parallel arrays of cylinders 106 (FIG. 1) to adjust the distance 109 (FIG. 1) between the arrays. One or both of the parallel arrays of cylinders 106 may be mounted on a moveable structure that may be moved to control the distance 2h 106 between the arrays. In this way, the structure may be tuned for input and/or output characteristics. In some implementations, output sensors may be used to determine one or more characteristics of the output of the system, which may be used by control circuitry and/or a processor to adjust the structure to provide the desired results.

It should be emphasized that the above-described embodiments of the present disclosure are merely possible examples of implementations set forth for a clear understanding of the principles of the disclosure. Many variations and modifications may be made to the above-described embodiment(s) without departing substantially from the spirit and principles of the disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

It should be noted that ratios, concentrations, amounts, and other numerical data may be expressed herein in a range format. It is to be understood that such a range format is used for convenience and brevity, and thus, should be interpreted in a flexible manner to include not only the numerical values explicitly recited as the limits of the range, but also to include all the individual numerical values or sub-ranges encompassed within that range as if each numerical value and sub-range is explicitly recited. To illustrate, a concentration range of “about 0.1% to about 5%” should be interpreted to include not only the explicitly recited concentration of about 0.1 wt % to about 5 wt %, but also include individual concentrations (e.g., 1%, 2%, 3%, and 4%) and the sub-ranges (e.g., 0.5%, 1.1%, 2.2%, 3.3%, and 4.4%) within the indicated range. The term “about” can include traditional rounding according to significant figures of numerical values. In addition, the phrase “about ‘x’ to ‘y’” includes “about ‘x’ to about ‘y’”. 

1. A system for generation of a higher harmonic of monochromatic electromagnetic radiation, comprising: a double array of sub wavelength dielectric cylinders, comprising: a first array of periodically spaced parallel cylinders, the periodic spacing corresponding to a resonant frequency; and a second array of periodically spaced parallel cylinders, the second array in parallel with the first array; and a source configured to illuminate the first array with a beam of monochromatic electromagnetic radiation at the resonant frequency to produce a higher harmonic of the monochromatic electromagnetic radiation.
 2. The system of claim 1, wherein the frequency of the higher harmonic of the monochromatic electromagnetic radiation is twice the resonant frequency.
 3. The system of claim 1, wherein the periodic spacing is about a wavelength that corresponds to the resonant frequency.
 4. The system of claim 1, wherein the second array is separated from the first array by about a wavelength that corresponds to the resonant frequency.
 5. The system of claim 1, wherein the cylinders are made of a nonlinear dielectric material.
 6. The system of claim 5, wherein the nonlinear dielectric material comprises a linear dielectric constant ∈_(c)>1 and a second order susceptibility χ_(c)<<1.
 7. A system for generation of a higher harmonic of monochromatic electromagnetic radiation, comprising: a double array of planar structures, comprising: a first periodically structured planar structure, the first structure having a periodic spacing corresponding to a resonant frequency; and a second periodically structured planar structure, the second structure in parallel with the first structure; and a source configured to illuminate the first structure with a beam of monochromatic electromagnetic radiation at the resonant frequency to produce a higher harmonic of the monochromatic electromagnetic radiation.
 8. The system of claim 7, wherein the second structure is separated from the first structure by about a wavelength that corresponds to the resonant frequency.
 9. The system of claim 7, wherein the second structure is separated from the first structure by less than a wavelength that corresponds to the resonant frequency.
 10. The system of claim 7, wherein the frequency of the higher harmonic of the monochromatic electromagnetic radiation is twice the resonant frequency.
 11. The system of claim 7, wherein a distance between the first and second periodically structured planar structures is adjustable.
 12. The system of claim 7, wherein the double array of planar structures comprises a double array of dielectric cylinders.
 13. A method for generation of a higher harmonic of monochromatic electromagnetic radiation, comprising: directing a beam of monochromatic electromagnetic radiation at a resonant frequency at a double array of planar structures, the double array of planar structures comprising: a first periodically structured planar structure, the first structure having a periodic spacing corresponding to a resonant frequency; and a second periodically structured planar structure, the second structure in parallel with the first structure; and generating a higher harmonic of the monochromatic electromagnetic radiation.
 14. The method of claim 13, further comprising adjusting a distance between the first and second periodically structured planar structures.
 15. The method of claim 14, further comprising sensing the generated monochromatic electromagnetic radiation, wherein the distance between the first and second periodically structured planar structures is adjusted based at least in part upon the sensed monochromatic electromagnetic radiation.
 16. The method of claim 13, further comprising focusing the beam of monochromatic electromagnetic radiation on the first periodically structured planar structure.
 17. The method of claim 13, wherein the double array of planar structures comprises a double array of dielectric cylinders.
 18. The method of claim 17, wherein the cylinders are made of a nonlinear dielectric material.
 19. The method of claim 13, wherein the frequency of the higher harmonic of the monochromatic electromagnetic radiation is twice the resonant frequency. 